/*==================================================
Project:       Targeting Social Programs
Authors:       Diether W. Beuermann
               Bridget Hoffmann        
               Marco Stampini 
               David L. Vargas
               Diego Vera-Cossio
----------------------------------------------------
Creation Date:    May 2023
Modification Date:   
Do-file version:    01
References:          
Output:             
==================================================*/

/* Regular PMT model using all SISBEN varibles and SISBEN income
 This will serve as baseline for comparison		*/


/*==================================================
            0: Program set up
==================================================*/
*Written on STATA 17
drop _all

/*==================================================
            1: load and transformations
==================================================*/

*----------- 1.0.1 Load CRRA welfare FUN
qui do "${dir6r}/ados/bstrap_pmt.ado" // load bstap FUN
qui do "${dir6r}/ados/bstrap_pmt_mh.ado" // load bstap FUN
qui do "${dir6r}/ados/baseline_pmt.ado" // load baseline FUN
qui do "${dir6r}/ados/expanded_pmt.ado" // load expanded FUN

*----------  1.1. Data prep:
use "${dir3r}/01_survey/survey_targeting_r1_long.dta", clear

* new interactions 
cap drop d_loss_illness_inc
cap drop d_reco_ilness_inc
gen d_loss_illness_inc = illness_incidence * d_loss_ilness
gen d_reco_ilness_inc = lag_illness_incidence * d_reco_ilness
lab var d_loss_illness_inc	"Loss: Illness Incidence"
lab var d_reco_ilness_inc	"Recovery: Illness Incidence"

cap drop d_loss_cut_remittance_inc
cap drop d_reco_cut_remittance_inc
gen d_loss_cut_remittance_inc = cut_remittance_incidence * d_loss_cut_remittance
gen d_reco_cut_remittance_inc = lag_cut_remittance_incidence * d_reco_cut_remittance
lab var d_loss_cut_remittance_inc "Loss: Remittances Cut Incidence"
lab var d_reco_cut_remittance_inc "Recovery: Remittances Cut Incidence"

* Globals with shocks
glo loss d_loss_ilness d_loss_cut_remittance death divorce bankrupt nat_shock
glo loss_rec      d_loss_ilness d_reco_ilness d_loss_cut_remittance d_reco_cut_remittance death divorce bankrupt nat_shock 
glo loss_inc      illness_incidence cut_remittance_incidence death divorce bankrupt nat_shock 
glo loss_rec_inc  d_loss_illness_inc d_reco_ilness_inc d_loss_cut_remittance_inc d_reco_cut_remittance_inc death divorce bankrupt nat_shock 

*Globals with updted survey variables
glo hh_char_srv urban age_feb20_srv prop_kids_srv n_members_srv elementary_srv highschool_srv tertiary_edu_srv children_srv living_couple
glo assets_srv own_washing_machine_srv own_tractor_srv own_moto_srv own_car_srv own_PC_srv own_fridge_srv
glo hh_srv  wall_finished_srv floor_finished_srv cookpower_connected_srv wc_flush_srv kitchen_srv electricity_srv running_water_srv trash_service_srv 
glo labour_srv formal_work_srv informal_work_srv

// Including labour status from SISBEN
glo pmt_1_srv  $hh_char_srv $assets_srv $hh_srv $labour_srv
** make sure all labels are OK 
lab var illness_incidence 			"Illness Incidence"
lab var cut_remittance_incidence 	"Remittances Cut Incidence"
lab var death						"Death"
lab var divorce					"Separation of Spouses"
lab var bankrupt					"Bankruptcy or Business Closure"
lab var nat_shock					"Natural Disaster or Fire"
lab var job_loss                   "Loss: Job lost at any point of the year"
lab var job_loss_formal            "Loss: Job lost formal at any point of the year"
lab var job_loss_informal          "Loss: Job lost informal at any point of the year"

*----------  1.1.2 Estimate PMT for survey data:


//============== Dynamic estimation 

// -------- OLS
local nreps=1000

//define poverty line in terms of Colombian Pesitos
local base=1.90*1515.1
// 2019 PPP factor https://data.worldbank.org/indicator/PA.NUS.PRVT.PP?locations=CO

*** Simple changes to changes
bstrap_pmt,  incvar(d_income_nt) variables(d_work_any_srv) reps(`nreps')  base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)
mat estmat0 = (1\1\1), r(estt)

*** Simple changes to changes with formal and informal 
bstrap_pmt,  incvar(d_income_nt) variables(d_formal_work_srv d_informal_work_srv) reps(`nreps')  base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)
mat amat = (2\2\2), r(estt)
mat estmat0 = estmat0 \ amat

*** Gain and losses // Main specification!
bstrap_pmt,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)
mat amat = (3\3\3), r(estt)
mat estmat0 = estmat0 \ amat

*** Gain and losses w/ formal and informal 
bstrap_pmt,  incvar(d_income_nt) variables(d_gain_formal_work_srv d_loss_formal_work_srv d_gain_informal_work_srv d_loss_informal_work_srv) reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)
mat amat = (4\4\4), r(estt)
mat estmat0 = estmat0 \ amat

*** Gain and losses at. def
bstrap_pmt,  incvar(d_income_nt) variables(d_gain_work_any_srv job_loss) reps(`nreps')  base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)
mat amat = (5\5\5), r(estt)
mat estmat0 = estmat0 \ amat

*** Gain and losses at. def w/ formal and informal 
bstrap_pmt,  incvar(d_income_nt) variables(d_gain_formal_work_srv job_loss_formal d_gain_informal_work_srv job_loss_informal) reps(`nreps')   base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)
mat amat = (6\6\6), r(estt)
mat estmat0 = estmat0 \ amat

*** Correcting for moral hazard

bstrap_pmt_mh,  incvar(d_income_nt) variables(d_gain_work_any_srv d_loss_work_any_srv) reps(`nreps')  mhazard(0.087) gainvars(d_gain_work_any_srv) lossvars(d_loss_work_any_srv) base(`base')  pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)
mat amat = (7\7\7), r(estt)
mat estmat0 = estmat0 \ amat

bstrap_pmt ,  incvar(d_income_nt)  variables(d_gain_work_any_srv d_loss_work_any_srv $loss) reps(`nreps') base(`base') pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv) 
mat amat = (8\8\8), r(estt)
mat estmat0 = estmat0 \ amat

*** using the cross-section PMT with baseline data but updating employment:

baseline_pmt, pmtvars(l_pp_inc_srv $pmt_1_srv)  predvars($assets_srv) negonly(no) reps(`nreps')  base(`base')
mat amat = (9\9\9), r(estt)
mat estmat0 = estmat0 \ amat

*** Main models

baseline_pmt, pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)  reps(`nreps')  base(`base') negonly(no)
mat estmat_static = (1\1\1), r(estt)
baseline_pmt, pmtvars(l_pp_inc_srv $pmt_1_srv)  predvars($labour_srv) negonly($assets_srv) reps(`nreps')  base(`base')
mat amat = (2\2\2), r(estt)
mat estmat_static = estmat_static \ amat
expanded_pmt,  pmtvars(l_pp_inc_srv $pmt_1_srv) predvars($pmt_1_srv)  reps(`nreps')  base(`base') 
mat amat = (3\3\3), r(estt)
mat estmat_static = estmat_static \ amat

/*==================================================
            4: Consolidate data
==================================================*/

** send to stata
cap frame create points
frame change points 

** append matrices

** dynamic points
// OLS
drop _all 
svmat estmat0, names(col)
gen type_m = "OLS"
g model_n = 3 + c1 if c1 != .
tempfile dfdyn
save `dfdyn', replace


// Static models

drop _all 
svmat estmat_static, names(col)
g model_n = c1 
gen type_m = ""
replace type_m =  "Baseline"  if model_n == 1
replace type_m =  "Updated"   if model_n == 2
replace type_m =  "Expanded"  if model_n == 3
append using `dfdyn'
tempfile dfdyn
save `dfdyn', replace


cap drop base _merge 
drop c1

foreach v in covered inclusion_err inclusion_ok exclusion_err exclusion_ok zero epov {
    replace `v' = `v' * 100
    replace `v'_uci=`v'_uci*100
    replace `v'_lci=`v'_lci*100
}

foreach v in covered hh_benefits inclusion_err inclusion_ok exclusion_err exclusion_ok crra crra_l crra_u zero epov covered_uci covered_lci inclusion_err_uci inclusion_err_lci inclusion_ok_uci inclusion_ok_lci exclusion_err_uci exclusion_err_lci exclusion_ok_uci exclusion_ok_lci crra_uci crra_lci crra_l_uci crra_l_lci crra_u_uci crra_u_lci hh_benefits_uci hh_benefits_lci budget_nat_uci budget_nat_lci zero_uci zero_lci epov_uci epov_lci {
    sum `v' if year == 2019 & model_n == 1
    replace `v' = r(mean) if year == 2019 & model_n != 1 
}

 // PPP values 

gen hh_benefits_lcu = hh_benefits
replace hh_benefits = hh_benefits / 1515.1 
replace hh_benefits_uci= hh_benefits_uci/1515.1
replace hh_benefits_lci= hh_benefits_lci/1515.1
// 2019 PPP factor https://data.worldbank.org/indicator/PA.NUS.PRVT.PP?locations=CO


// store data for later use 
save "${dir3r}/03_auxiliar/Boostrap_estimates.dta" , replace
